Analytic Lyapunov exponents in a classical nonlinear field equation 
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It is shown that the nonhnear wave equation di (p — d^(p — 
tJ-odx{dxip)'^ ~ 0, which is the continuum limit of the Fermi- 
Pasta-Ulam (FPU) /3 model, has a positive Lyapunov expo- 
nent Ai, whose analytic energy dependence is given. The re- 
sult (a first example for field equations) is achieved by evalu- 
ating the lattice-spacing dependence of Ai for the FPU model 
within the framework of a Riemannian description of Hamil- 
tonian chaos. We also discuss a difficulty of the statistical 
mechanical treatment of this classical field system, which is 
absent in the dynamical description. 

PACS: 05.45. +b; 03.40.-t; 03.40.Kf 



The numerical study that E.Fermi, J. Pasta and S.Ulam 
made at Los Alamos, on a chain of anharmonically cou- 
pled oscillators, represents a milestone in the develop- 
ment of the modern theory of nonlinear dynamics [Q. 
The surprising outcomes of this first computer experi- 
ment in the history of physics have been very seminal. 
Among the other attempts to explain the apparent lack 
of thermalization of the energy initially stored in one of 
the linear modes of the chain, Zabusky and Kruskal [^|J|] 
remarked that a continuum limit version of the Fermi- 
Pasta-Ulam (FPU) model leads to the Korteweg-de-Vries 
(KdV) equation, when the a model is considered, and to 
a modified KdV equation when the /3 model is consid- 
ered. These are nonlinear, integrable partial differential 
equations where a special class of solitary waves - that 
Zabusky and Kruskal called solitons - exists and can to 
some extent explain the FPU recurrences Q . 

In the present work we tackle the FPU-/? model, de- 
scribed by the Hamiltonian 
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and we consider another legitimate continuum limit |^] 
of this system, described by the Hamiltonian 



H{tt, 



dx 



i-7r(.)^ + |(V0)^-f^(V# 



(2) 



and leading to the field equation (setting rj = v = 1) 



Also this continuum limit of the FPU model might ap- 
pear integrable on the basis of the following reasoning 
Let us consider the Legendre transform: x — dxff'j 
dt(t> and x — d^^ip, t = i9rV' with the relation 
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where l/D = d^ipd^tp — (d^^-ip)"^ . 



Eq.(||) we can linearize it to 



Substituting into 
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which can be solved by variable separation by setting 
-4, = F{x)G{t), where G" ± c^G = 0, and F" ± ^{l + 
?i11qX^)F = 0. In principle, by inverting this Legendre 
transform, Eq.(^) could be analytically solved. However, 
as we shall discuss, the method sketched above can only 
lead at most to local and not to global invertibility, both 
in time and space. Actually we are going to show that the 
field equation (|^) for the FPU (i model in the continuum 
limit is chaotic (the Lyapunov characteristic exponent 
of the solutions is positive). The pattern of the Lya- 
punov exponent that we shall derive is shown in Fig.l. 
The cross-over in the energy dependence of the largest 
Lyapunov exponent has been attributed in Refs. to a 
(smooth) transition between weak and strong chaos 
The persistence of such a cross-over also in the continuum 
limit is a remarkable fact. Loosely speaking, at high en- 
ergy, in the strongly chaotic regime, the spatio-temporal 
behavior of the field (t>{x, t) should look like a fully de- 
veloped "turbulent" field, whereas this should not be the 
case at low energy, where a weakly "turbulent" dynam- 
ics would set in. Our calculation shows that, no matter 
whether weak or strong, chaos is present in the system at 
any energy, and since the system is Hamiltonian, we can 
believe that a standard statistical mechanical description 
is allowed. This raises an important problem that we now 
discuss. 

For the lattice 

system corresponding to (^ the standard canonical par- 
tition function Z = j Y\f=i dTTid4>i exp [—f3H{{'Ki}, {^i})] 
in the limit of lattice spacing a ^ becomes 



Z = 



(5) 



It is interesting at this point to raise a problem that 
comes out when some expected dynamical properties of 



1 



the field (j){x, t) are compared to their statistical counter- 
parts worked out by means of Z. By means of an orthog- 
onal change of coordinates, the Hamiltonian (|l|) is rewrit- 
ten in the form § H{Pk,Qk) =Y.kEk= EJll^l + 
^kQl)+f^J2kiM,k3 '^(^' ^1' ^2, h)QkQkiQk^Qk3] where 
the summations run from k — 1 io k — k„iax for 
a lattice, and from ko to k„iax = oo in the contin- 
uum limit. In the new variables the partition function 
becomes Z = ^ WlZT dPkdQke^^[-PH{{Pk},{Qk})] 
which, in the limit kmax ^ oo, gives Z = 
J VP{k)VQ{k) exp \-f3 dkH{P{k), Q(fc))j , 

where fcg = ^ > for a finite support. The generalized 
equipartition theorem states that {Q^dH/dQu) = 
const independently of fc, and from Ref. | [Tc[ | we know that 
for the FPU /3 model such an average is (QkdH/dQk) — 
k^Q^{l + a), where a is a constant, and therefore ~ 
k^'^{l + a)^^, which, being independent of k^ax, still 
holds true in the continuum limit, i.e ~ fc~^(l -f- 

a)"^. Whence an ultraviolet catastrophe, i.e. the diver- 
gence of average total energy when k^ax oo also on a 
finite support [0,L]. However such a difficulty is absent 
in the dynamical equation (^, which, being derived from 
a Hamiltonian, must keep the energy constant, therefore 
finite if it is finite at time t = 0. 

The point is that statistical mechanics and dynamics 
yield very different large-k spectral properties. In fact, 
already on the basis of the purely dimensional analysis 
discussed in Refs. [H, we can easily find that in order to 
bound the total energy (on a finite support [0,L]) to a 
finite value, the ultraviolet asymptotic power spectrum 
of (f){x) must be bounded above by ~ which 

means that the high spatial-frequency modes cannot have 
the same energy of the low spatial-frequency modes. 

The rough estimate of the ultraviolet spectrum of 
the continuum model can be improved as follows. Any 
smooth initial condition (j)[x^ 0) S C°°(M) of the equation 
of motion (|^) has to remain smooth at any further time 
t, i.e. (I){x,t) e C°°(R). Now, the analytic continuation 
^{z, t) of (p{x, t) will be of class on a strip whose width 
i5 will be determined at any given time t by the distance of 
the closest singularity to the real axis. If z^^ ~ xf^+iy^ are 
the locations of the singularities of ^{z,t), to the lead- 
ing order one finds $(/c) - /c"" 'Zn Cne'''< for 
/c — > oo, where, apart from the power-law prefactor 
the sum is carried over the Imz > half plane and c„ 
are suitable coefficients. The factors e*'^^" = e*'^^'«e~'°^" 
efficiently single out the closest singularity to the real 
axis to give |$(fc)p < Ae~^'^ at k ^ oo which provides 
a better ultraviolet estimate of |0(A:)p stemming from 
the constraint of differentiability of (j>{x). Needless to 
say, this ultraviolet exponential fall off of |(?I>(fc)p guar- 
antees the finiteness of the total energy for any finite 
support [0,L]. It is not unreasonable to think that a 
relationship might exist between the minimum width of 
the spatio-temporal cell where the Legendre transform 



of Eq.(||) - sketched at the beginning of this paper - is 
invertible, that is where Eq.(^) is locally integrable, and 
the width of the {k,uj{k)) cell where the power spectrum 
|(/)(A:)P exponentially falls off. There is apparently no 
way of reproducing the ultraviolet exponential decay of 
|(/)(A:)P - which is naturally brought about by the dy- 
namical equation ^ - within the statistical mechanical 
framework. The field 4>{x, t) can be chaotic in space and 
time only down to some small scale which naturally pro- 
vides a cutoff that, inserted into Z, would remove the 
mentioned divergences. Let us mention that the physical 
meaning of a rigorous treatment of another classical non- 
linear field equation (complex Ginzburg-Landau) [ p^ is 
coherent with the theoretical scenario depicted above. 

To come now to our main result, chaoticity in the con- 
tinuum limit, we calculate the Lyapunov exponent by 
extending to such a limit a method to tackle Hamilto- 
nian chaos that has already given excellent analytic pre- 
dictions of the largest Lyapounov exponent in the ther- 
modynamic limit of the FPU /3-model on a lattice. 
This method exploits the mathematical identification of 
an Hamiltonian flow with the geodesic flow on a suit- 
able Riemannian "mechanical manifold" consisting of an 
enlarged configuration spacetime ({g° = t,q^, . . . , q'^} 
plus one real coordinate (7^+^), whose arc- length ds^ = 
-2V{q){dq°)^ + a^jdq^dq^ + 2dq°dq'^+^ defines the so- 
called Eisenhart metric |l^]; V is the potential. In 
the geometrical framework, the (in)stability of the tra- 
jectories is the (in)stability of the geodesies, and it is 
completely determined by the curvature properties of 
the underlying manifold according to the Jacobi equa- 
tion for the geodesic deviation. This equation, writ- 
ten for Eisenhart metric, entails the usual tangent dy- 
namics equation -I- (d'^V/dqidq^)^^ = 0, which is used 
to measure Lyapunov exponents in standard Hamilto- 
nian systems. Having recognized its geometric origin, it 
can be transformed, under certain hypotheses [ p^ , into 
an effective scalar stability equation that, independently 
of the knowledge of dynamical trajectories, provides a 
measure of their average degree of instability. This ef- 
fective stability equation is in the form of a stochas- 
tic oscillator equation -ip + [kq + ai^ri{t)]ip — [p^ , 
where = {Ki^}/N, al = {{Kr - {Kr))^)/N , and 
Kji = AV = J2iLi^'^^/^1i Ricci curvature of 

the mechanical manifold, computed with gE] ri{t) is a 
gaussian (5-correlated random process of unit variance. 
The exponential growth rate A of {tp'^ + "0^) is computed 
exactly to provide the following estimate of the largest 
Lyapunov exponent: 



A = 
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where t = 7rY/7to/(2[Ko('«o + o'k)]^ -|- ttCTk) is a charac- 
teristic time scale worked out on the basis of geometrical 
arguments. In the limit (Tk/kq ^ 1 one finds A oc cr^ [p5[. 
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In this geometric picture chaos is mainly originated by 
the parametric instabihty activated by the fluctuating 
curvature "felt" by the geodesies. 

The quantities kq and cr^ can be exactly computed for 
the FPU-/? model as microcanonical averages by taking 
advantage of the analytically known canonical partition 
function and by using the conversion formulas relating 
canonical and microcanonical averages. The final ana- 
lytic expressions are given in Ref. , and among them 
we report those needed here 
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where A{d) 



D_ 



2(e) 



D 



-1/2 



with -D_3/2 and -D_i/2 parabolic 



cylinder functions, the parameter is a function of the 
inverse temperature /3 through 9 = and e{6) is 

the energy per degree of freedom of the system. The mi- 
crocanonical average of curvature fluctuations involves a 
correction term to their canonical average which involves 
the specific heat Cv{9) computed, as usual, as the second 
derivative of the free energy. These quantities enter the 
formula (||) to give the analytic values of Lyapunov expo- 
nents at different energies. Now, it is useful to work out 
simple analytic expressions of Ai (e) at low and high en- 
ergy densities. In order to obtain them, notice that from 
Eq.(^ one gets /ie = f{9), a function of 9; this func- 
tion is invertible, because in the absence of phase transi- 
tions e(T) ( T is the temperature ) is always one-to-one. 
Therefore we can write 9 — /~^(/ie) and consequently 
Ko and tr^ can be expressed as functions of /ie. Making 
use of the asymptotic expressions of the parabolic cylin- 
der functions, this inversion is easily obtained in the two 
opposite limits e ^ and e — > cxd. Also simple analytic 
expressions for kq, and r can be obtained in these two 
limits. At e ^ 0, and 9 ^ 1, we get 
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(10) 
(11) 



Let us remark that there is some arbitrariness in the 
derivation of r in Ref. |l^; however from Eq.(^ we see 
that the only relevant constraint on r is that ct^t — > 
for e ^ 0, so that A - and Ai with it - can be expanded 
in series of the small quantity d^r yielding 



(13) 



which is in strikingly good agreement with the low energy 
"experimental" results for Ai. Similar developments can 
be worked out in the limit e oo (and so 9 1) giving 
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(14) 



by computing again T(e) we find T(e) ~ c(/ie) yield- 
ing 

c'(Hi/4 



Ai(e) 



(15) 



where c and c' are constants. Also this scaling law is in 
perfect agreement with numerical results. The compari- 
son between the full analytic prediction of Ai(e), obtained 
by substituting Eqs.(||) and (||) into Eqs.(^), shows a per- 
fect agreement with the outcomes of a standard numer- 
ical computation of the largest Lyapunov exponent at 
different values of e. Let us remark that the analytic pre- 
diction of Ai(e) is obtained in the N oo limit, which is 
implicit in the computation of the microcanonical aver- 
ages of the Ricci curvature and of its fluctuations. Now, 
as the analytic expression in (^) is in excellent agreement 
with the numerical data obtained for discrete (particle) 
systems at different values of N, there is no reason to 
doubt that such an agreement will hold true also when 
N is varied while keeping the total length L of the sys- 
tem fixed, i.e. in the continuum limit. The non-trivial 
difference between this limit and thermodynamic limit in 
Ref. [|l^ is that N ^ oo while the energy remains finite 
in the former case, whereas both N and energy diverge 
in the latter case. Thus let us now consider the Hamilto- 
nian ^ discretized on a lattice of length L and spacing 
a (with L = Na) 

H{{n,},{c^,})^ 

N 
i=l 

where the dimensional constants i^, rj have been intro- 
duced. We have thus obtained a discrete (lattice) system, 
similar to that described by the Hamiltonian (|l|) ; in order 
to make a direct comparison between the two Hamiltoni- 
ans of Eqs.(p^ and (^, we must now absorbe the dimen- 
sional parameter a into the constants of the Hamiltonian. 
By means of the following rescaling of the coordinates 0i 
and of the time t: ~ 4>i/a^ with a — \fajr\^ t' — t/^ 
and 7 = \fr\Jv^ and denoting by \C,i\ the momenta con- 
jugated to K of Eq.(|l6|) is put in the form 
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which, with the explicit computation of r gives 



(17) 
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which is just the form of H in Eq.(|l|), so that we can ap- 
ply the geometric treatment of dynamical chaos also to 
the systems described by this family of Hamiltonians Ha ■ 
The precise meaning of this comparison is the following. 
The continuum Hamiltonian (Q) has H of Eq. ( ]l6| ) as 
its discretized version, where the parameters 77, v, /iq are 
equal to the values of the continuum model. After the 
recasting of Eq. (|l^) into the form (|l^) we can easily de- 
rive the value of the coupling constant /j, of the discrete 
FPU model ffl) so that a correspondence can be made be- 



tween the discrete model and the lattice version (17) of 



the continuum model for any value of the lattice spacing 
a. Evidently the relation between the coupling constants 
is /i = ^Q/arf' , and this tells that in order to represent the 
solutions of a set of systems (|l^) by means of those of a 
set of discrete systems (|l|) with a finer and finer sampling 
of the spatial support of length L, i.e. letting a ^ 0, the 
coupling constant n must be larger and larger. To the 
lattice-discretized system (|l^) we can now apply all the 
above equations for the geometric description of chaos, 
provided that we replace /i with /io/ary^. In passing to 
the continuum limit we have to let a — > while the total 
energy of the system Etot is kept fixed. Hence the energy 
density e has to diminish as Etot/N^ but N — L/a and 
then e = Etot ■ a/ L. Hence we obtain 

Mo -Etot a 



arj'^L 



cost 



(18) 



From Eq.QlSD it then follows that the equations for kq and 
(Tfj (||)-(^ are stable with respect to the limit to the con- 
tinuum and so is the largest Lyapunov exponent. There- 
fore the field equation (^ is chaotic and its Lyapunov ex- 
ponent has the pattern shown by Fig.l. This constitutes 
a first example of analytic calculation of Ai(e) for a field 
equation. Implicit in our computation is the possibility of 
approximating, as accurately as needed, a partial differ- 
ential equation by means of a finite (truncated) system of 
ordinary differential equations, for which Lyapunov expo- 
nents are well defined. In conclusion we have shown that 
the continuum limit of the Fermi-Pasta-Ulam /3 model is 
chaotic. Wc have discussed the interesting problematics 
raised by this result on the statistical description of clas- 
sical field theory models that might be relevant also to 
the Wick-rotated quantum field theories. 
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FIG. 1. Analytic Lyapunov exponent for the continuum 
limit as a function of the physical parameters. 
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